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Abstract Solar tomography has progressed rapidly in recent years thanks to 
the development of robust algorithms and the availability of more powerful 
computers. It can today provide crucial insights in solving issues related to the 
line-of-sight integration present in the data of solar imagers and coronagraphs. 
However, there remain challenges such as the increase of the available volume 
of data, the handling of the temporal evolution of the observed structures, 
and the heterogeneity of the data in multi-spacecraft studies. We present a 
generic software package that can perform fast tomographic inversions that 
scales linearly with the number of measurements, linearly with the length of 
the reconstruction cube (and not the number of voxels) and linearly with the 
number of cores and can use data from different sources and with a variety of 
physical models: TomograPy (http://nbarbey.github.com/TomograPy/), an open- 
source software freely available on the Python Package Index. For performance, 
TomograPy uses a parallelized-projection algorithm. It relies on the World Coor- 
dinate System standard to manage various data sources. A variety of inversion 
algorithms are provided to perform the tomographic-map estimation. A test 
suite is provided along with the code to ensure software quality. Since it makes 
use of the Siddon algorithm it is restricted to rectangular parallelepiped voxels 
but the spherical geometry of the corona can be handled through proper use 
of priors. We describe the main features of the code and show three practical 
examples of multi-spacecraft tomographic inversions using STEREO /EUVI and 
STEREO/COR1 data. Static and smoothly varying temporal evolution models 
are presented. 
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1. Introduction 
1.1. Motivation 

Except for the rare case when in- situ exploration is practicable, the properties 
of astronomical objects are deduced from the analysis of the properties of light 
only. Most astrophysical measurements are therefore affected by the problem of 
line-of-sight (LOS) integration, i.e. the modification of the signal of interest by 
background and foreground emission and absorption. This problem is one of the 
major sources of uncertainties in the diagnostics of the solar plasma. 

Integration along the LOS tends to confuse structures to the point that 
measurements crucial to the understanding of coronal physics are difficult to 
interpret. The controversy about the nature of polar plumes is one example. 
Polar plumes are observed at visible and UV wavelengths extending quasi- 
radially over the solar poles. Their appearance in photographic records led to the 
classical view of plumes as being pseudo-cylindrical structures denser than the 
surrounding corona. However such linear features in the images can also result 
from chance alignments of fainter structures distributed along a network pattern 
and integrated along the LOS. Both types of plumes have been supported by 
different authors, and it is possible that the two types coexist. See, e.g., |Gabriel| 



et al. (20091 for a detailed discussion. Since the two proposed types of plumes 
have nearly identical properties in remote-sensing data, the true nature of these 
objects remains a subject of controversy. 

One can also cite the problem of background estimation in the coronal-loop 
debate. As building blocks of the solar corona, loops have been extensively 
studied. However, we still cannot answer fundamental questions such as what 
processes are responsible for their formation or for their heating. One of the fac- 
tors explaining this state of facts is that the determination of physical parameters 
such as density and temperature within the loops is rendered difficult by LOS 



superimposition. For example, Terzo and Reale (2010) have shown that differ- 



ent estimations of the loop-background radiation lead to different temperature 
profiles and different conclusions regarding the loop cooling. 

Different strategies have been devised over the years to overcome the lim- 
itations imposed upon remote-sensing data by LOS integration. One obvious 
approach is to select an observation time when the solar corona presents a 
simple geometry for which it is possible to estimate the contribution of the 
various regions of the LOS. When observing polar plumes, for example, data are 
acquired preferentially at solar minimum when the polar holes are well developed 
so that the contribution of streamers to the foreground and background is min- 
imum. However, even if these conditions are met, it is likely that several plumes 
or plumes and inter-plumes are superimposed along the LOS, thus confusing 
the interpretation. Therefore, favorable observing conditions are generally not 
sufficient to exclude possible LOS ambiguity. 

If a simple coronal configuration cannot be assumed, or more generally if 
superimpositions cannot be ruled out , one has to devise means of analyzing the 
LOS content. Spectroscopic techniques such as the Differential Emission Measure 
(DEM) can be used to estimate the quantity of emitting plasma along the LOS 
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as a function of temperature. If this approach is able to detect the presence 
of regions of different temperatures along the LOS, it does not say how the 
temperatures are distributed spatially. A single multithermal volume can have 
the same DEM signature as the superimposition of several large-scale isothermal 
structures. 

Line-of-sight ambiguities can be alleviated, at least partially, if one can make 
several simultaneous observations from different locations. The twin spacecraft 



of the STEREO mission (Kaiser et al., 20081 were designed to achieve this 



The two vantage points that they offer provide precious information on the 
LOS content. In some cases, especially with a high-contrast object having well- 
defined boundaries such as coronal loops, direct stereoscopic reconstructions can 
be performed. Such reconstructions can, for example, be used to assess the qual- 



ity of the background estimation used in loop studies (e.g. Aschwanden et al 



2008). However for more diffuse objects not presenting sharp boundaries, or for 
which the visible boundaries can be LOS-integration artifacts, such as streamers, 
plumes, etc., direct, stereoscopic reconstruction is not reliable. However, if two 
viewpoints or more are available, tomography is a possible approach to inverting 
the LOS integration. 

1.2. Solar Tomography 

The term tomography encompasses a wide range of techniques aimed at imaging 
the internal structure of objects. Tomographic techniques are used in many 
areas of scientific research such as medicine, geophysics, materials science, and 
astrophysics. In the particular case of solar tomography, images recording the 
line-of-sight integration of coronal emission and taken from different viewpoints 
are used to estimate local physical quantities such as the electron number density 
or temperature. This is achieved using computed tomographic reconstruction 
techniques identical to the ones used in medical computer tomography. Math- 
ematically, it is an inversion of the line-of-sight integration. This method is 
sometimes called Solar Rotational Tomography (SRT) as it generally relies upon 
the solar rotation to simulate data acquisition from different viewpoints. 

However, there is significant differences between SRT and medical computer 
tomography that renders the problem more difficult to solve in the case of SRT. 
First, medical imaging scenarios benefit from high signal to noise ratio (SNR) 
measurements and have much higher measurement density than what is currently 
available with solar observatories. But more importantly, medical imaging has 
much less restriction on the number of point-of-view. Indeed, most of the time, 
SRT is restricted to one instantaneous point-of-view, the only exception being the 
STEREO twin spacecrafts. Another important difference with medical imaging 
is the presence of an opaque sphere in the middle of the region of interest: the 
photosphere. A similar issue is the presence of occulter in coronograph instru- 
ments, which restricts spatial information available in the data. This renders the 
problem much more ill-posed in SRT than in medical imaging tomography. 

Progress in solar tomography comes from the availability of new data, new 
physical models, new inversion algorithms, and more powerful computers. Solar 



tomography can be traced back to Van de Hulst (19501 who presented a one- 



dimensional inversion of white-light data by fitting an analytic corona model. 
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Another seminal paper of SRT is Altschuler and Perry (19721. It introduced 
computer aided numerical estimation of the three-dimensional electron density of 
the corona using data from the K-coronameter of the High Altitude Observatory. 
The method reduces to a least-square estimation of the coefficients of Legendre 
polynomials. 

Later, Davila (19941 investigated, through simulations, the possibility of per- 
forming full three-dimensional (3D) solar tomography of the corona. That article 
assumed that more than one spacecraft would take data (up to nine actually) 
and thus did not require data to be taken at different times. Emission-map 
estimation was performed using the algebraic reconstruction technique (ART) 
which is an iterative gradient method used to invert the linear tomographic 
model. Estimated maps where reduced to 20 x 20 grids. 

With the availability of the Solar and Heliospheric Observatory (SOHO) 
data came the first three-dimensional maps of the corona both in white light 
[e.g. Frazin and Janzenj 2002 1 with the Large Angle Spectrometric Coronagraph 



(LASCO: Brueckner et al. 19951 and in the ultraviolet (e.g. Panasyuk 1999) 



using the UltraViolet Coronagraph Spectrometer (UVCS: |Kohl et al. 1995). An 
algorithm-oriented article by Frazin (2000) introduced a penalized likelihood 
approach minimized using an iterative solver (conjugate-gradient) to allow noise 
mitigation through proper regularization and modeling of the outliers. 

Several generalizations were then developed. For example, jWiegelmann and] 



Inhester (2003) introduced a method for the joint estimation of the electron 



density and the magnetic field while Frazin, Vasquez, and Kamalabadi (2009) 
proposed a reconstruction of the local DEM from EUV images. 

The temporal evolution of coronal structures during data acquisition is one 
of the main issues in coronal tomography. Two different directions have been 



investigated to address this issue: either assume slow evolution (e.g. Butala et al. 



20101 or further restricting the possible evolution to specific structures (e.g. 



Barbey et al. 2008). 



1.3. Outline 



In this article we describe TomograPy, an open-source software package imple- 
menting the main desirable features in a generic solar tomography code (see, 
e.g., Frazin and Kamalabadi 2005 for a review): capability to use both EUV and 



white-light data to estimate the local electron density and temperature, modeling 
the temporal evolution of structures during data acquisition, and performing 
rotational tomography with multiple spacecraft, i.e. with STEREO data. 

In Section 2.1 we formulate mathematically the tomographic inversion prob- 
lem and introduce the notations used in the description of the code given in Sec- 
tion[3] Section [4] gives an overview of the numerical performance of the code while 
Section [5] describes three practical examples of tomographic reconstructions. 
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2. Tomographic Inversion 

2.1. Linear Inverse Problem 



The problem to invert can be expressed by Equation ([TJ where y is the data, A 
is called the physical model, x is the object map to estimate, and n is an additive 
noise (which we will assume Gaussian, independent and identically-distributed 
(iid)) 

y = Ax + n (1) 

The physical model represents all of the transformations that link the quantity 
to estimate (e.g. the local emissivity, electron temperature, electron density, etc.) 
to the data. It always includes the line-of-sight integration and the model of 
temporal evolution, even if this later is implicitly static. The physical model 
may also include the formation process of the observed lines if one wants to 
estimate a physical quantity such as the local electron density or temperature 
instead of the local emissivity. We restrict ourselves to cases where the data are 
a linear function of the unknown quantities to determine. It is worth noting, 
however, that tomographic inversion can be done even if the physical model is 
non-linear, although accompanied by an important increase in complexity. An ex- 
emple of non-linear tomography application is electrical capacitance tomography 



Soleimani and Lionheart (2005) which is intrinsically non linear. Using Monte 
Carlo Markov Chain (MCMC) methods, it would be feasible to fit non-linear 
models of Coronal Mass Ejections as the one provided in |Thernisien, Vourlidas,| 



and Howard (2009). 



2.2. Bayesian View on Linear Inversion 



In the Bayesian paradigm, a probability density function (PDF) is associated 
with each variable. To invert the class of problems described by Equation (JT|), 
one needs to know the statistical properties of the noise [n], which gives the 
likelihood. One also needs to define a prior on the unknowns [x]: it gives the 
PDF of x knowing the data [y], which is called the posterior on x. This is done 
through Bayes' rule given in Equation ^ 

where A4 regroups all of the assumptions on the model. In this article, the PDF 
of x given y is noted f(x\y). 

In the case of a Gaussian multivariate likelihood and prior, the posterior is 
also Gaussian, and thus fully determined by its mean and covariance matrix. 
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This is summed up in Equation ([3| 

f(x\M) ~ N{V,al(B T B)- 1 ) 

f(y\x,M) ~ M{Ax,a 2 n I) 

f(x\y,M) ~ 7V(x,±) 

x = i±A T y (3) 

E = (A T A + \B T B)~ 1 

A = "f 

where Af(ji, E) is a multivariate Gaussian of mean and covariance S. In 
Equation ([3]) we assumed an independent Gaussian noise of variance a\, so 
that the covariance of the likelihood is cr 2 J where i" is the identity matrix. We 
defined a zero-mean prior [/(a;|Al)] with B being a prior model. B can be, for 
instance, a finite-difference operator. A zero-mean prior combined with a finite- 
difference operator means that the finite difference of the map tends to be close to 
zero. In other words, this prior would favor smoother solutions over non-smooth 
solutions close to the data. This is a sensible choice for electron density at the 
scales considered. a n and consequently A are free parameters. It is possible to 
assign a PDF to A in order to estimate this parameter in an unsupervised way 
but it results generally in very resource-consuming algorithms. In this article, 
we will restrict ourselves to a fixed A in a supervised way. 

Characterizing the solution in a Bayesian way requires the estimation of both 
S and x. x gives the most probable solution and X gives information about the 
uncertainties in the unknowns. However, in most practical cases, the covariance 
matrix [X] is too large to be stored in memory, and one only keeps x. In this 
case, a full matrix inversion is not required, and one can estimate x much faster 
using iterative schemes such as the conjugate gradient method. 

Since x is the maximum a posteriori (MAP) of the problem, it is also the 
minimum of the co-log-likelihood as written in Equation Q 

x = argmin { J(x)} 

X 

= argmin {- log [f(x\y,M)}} 

X 

= argmin {- log [f(y\x,M)f(x\(M))}} (4) 
.r 

= argmin |||y — Aa;|| 2 + A ||Ba;|| 2 | 
= (A T A + XB T By 1 A T y 

The term \\y — Ax\\ 2 is a simple least-squares term. It defines the closeness 
to the data. The second term [|jSa;|j] is a regularization term that prevents the 
estimate from being noisy. In terms of matrix inversion, A T A is ill-conditioned 
and B T B is added in order to have a better-conditioned matrix. To find x, 
iterative-gradient methods need only the definition of the criterion J(x) and its 



SOLA: siddon.tex; 18 January 2013; 20:30; p. 6 



TomograPy: Solar Tomography 



gradient V x J. Gradient methods can be order of magnitudes faster than the full 
inversion of the matrix, especially when A and B are sparse or when the problem 
has been properly preconditioned. 



3. Main Features of TomograPy 

3.1. Fast Parallelized Projector 



ica 



TomograPy is a Python (Van Rossum and Centrum voor Wiskunde en Informat- 



1995 ) package build around a C implementation of the Siddon algorithm 
(Sid don, 1985[ ) of line-of-sight integration, and is thus restricted to rectangular 
parallelepiped voxels. This C projector has been parallelized using OpenMP 



(Dagum and Menon, 2002). Numpy (Oliphant, 2006) is a requirement as well 



as PyFITS (Barrett and Bridgman, 



ally, one can use SciPy (Jones, Oliphant, and Peterson, 2001 



1999|) to handle FITS data files. Option- 
sparse matrix 



optimization routines to perform fast linear inversions. The algorithm has been 
carefully optimized using meta-programming techniques to avoid if statements 
and function pointers in the inner loop. This has been done using templates of 
C code, and replacing key values in the source template to generate variations 
in the C sources for various application (e.g.: float and double values, projection 
and backprojection, presence of an obstacle or not). Here the word template is 
not to be confused with C++ templates but is more closy related to the notion of 
web template. The same idea is used in Numpy itself and allow more flexibility 
to pure C code. 

The projection algorithm provided with TomograPy can be used with a variety 
of estimation methods as long as they rely on the linear-operator interface. It 
allows for fast testing of various optimization strategies. Results presented in this 
article will exclusively use conjugate gradient schemes, but TomograPy provides 
other options. The key requirement is for the algorithm to rely on matrix- vector 
operations. 

See Section [4] for an analysis of the performance and scaling of this implemen- 
tation. The TomograPy projector is well tested and provided with a test suite, 
which fully covers this part of the code. 



3.2. Instrument Independence 



TomograPy takes as input FITS files (Flexible Image Transport System: Wells, 
Greisen, and Harten 1981 1 containing fully calibrated images expressed in units 



consistent with the physical model chosen by the user. TomograPy internally uses 



the World Coordinate System (WCS) ( |Calabretta and Greisen, 2002 Greisen 
and Calabretta, 2002 ) standard keywords to determine the position of the ob- 



server and to define the projector from the data and the desired format of the 
object map. TomograPy will therefore accept any data compliant with the WCS 
standard. As the data of most current instruments are already provided as WCS- 
compliant FITS files, all that is required is to store a set of calibrated files in 
a directory that TomograPy will be pointed to. For data that does not conform 
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to WCS, it is straightforward to write a small wrapper that will handle the 
instrument-specific metadata and convert them into the corresponding WCS 
keywords. 

TomograPy allows inversions using data from multiple spacecraft, for example 
with STEREO- A and B and SOHO. The data from the different instruments 
nonetheless need to be consistent, i.e. to record the same physical quantity. 

3.3. Physical models 

As stated in section |2.1| any linear model can be inverted using the same 
framework. TomograPy provides with the possibility to perform inversions with 
the several models described in this section. The TomograPy projector already 
discussed in Section |3.1| is a building block for all of the models described here. 
We will first describe models of temporal evolution and then models of line 
emission. It is possible to combine these models to perform, for example, multi- 
spacecraft, smooth, temporal rotational tomography. In the future, it will be 
possible to combine the models presented here with models not yet available in 



TomograPy such as the Differential Emission Measure model (Frazin, Vasquez, 



2003). 



and Kamalabadi, 2009 1 or magnetic-field models ( Wiegelmann and Inhester, 



3.3.1. Single Spacecraft Static Tomography 

This is the simplest case. In the next two sections we will consider only the line- 
of-sight inversion without assumptions on the line formation process. In this case, 
static rotational tomography of the solar corona can easily be formulated as in 
Equation once discretized. In this article, we will assume that the object-map 
cube has been discretized using contiguous rectangular parallelepiped voxels of 
identical shape; it is a requirement of the Siddon algorithm. Since the intensity 
on one detector results from the integration along the line of sight of the emission 
in the observed object, it can be expressed as in Equation ([5]): 

yj = ^ p iJ e i + n i ( 5 ) 

i 

where yj is the intensity of the detector j , ej is the emission in the voxel j , and 
Pi j is the length of the segment of the line of sight i that corresponds to the 
voxel j, and rij is the noise observed on detector j . Reformulating Equation ^ 
in terms of vectors and matrices, we obtain Equation ([6| where the t index refers 
to the time at which the data have been taken. 

Vt = Pt.e t + n t (6) 

where P, the TomograPy projector, is the most basic block for building physical 
models. Note that we indexed both axes of the detector using one index as well 
as the voxels of the object map. We can do the same on a time index, since we 
assume that there is no temporal evolution (e does not vary with t). Regrouping 
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all of the data at all considered instants [t], results in Equation ^ which is 
similar to Equation ([!]). 



(7) 



This model can be used to estimate emission maps from EUV data using 
formula Q. We can use a smoothness prior to avoid having too much noise 
in the maps. This is done using a finite-difference operator along each axis of 
the maps for the B matrix. If we want to account for the lower signal-to-noise 
ratio that we typically have in solar rotational tomography, we can have finite- 
difference operators weighted by the altitude of the considered voxels. Finally, 
we have the following equation 

e S RT = argmin {\\y - Pef + A \\DRe\\ 2 } (8) 

where D is the finite-difference operator and R is a diagonal operator with the 
height of the voxels on the diagonal. The use of a smoothness prior increasing 
linearly with height, allows for the maps not to be affected by the difference 
between spherical grids and Cartesian grids. Spherical grids have bigger voxels 
at high altitudes increasing the SNR per voxel with height. This is not the case 
for Cartesian grids but it is compensated by the use of a height-dependent prior. 

3.3.2. Multiple- Spacecraft Static Tomography 

If now we want to use data from multiple spacecraft, with the static assumption 
we can use Equation (J9j) (assuming two spacecraft A and B without loss of 
generality) . 

Va = 9APAe + n A 

Vb = gBPBe + n B (9) 

Each of Equations (|9| are derived from ^ but both spacecraft can have 
different gains (g^ and gs) at the wavelength considered. In this model, it is 
not possible to assume a different behavior of the filters as a function of the 
wavelength since e needs to correspond to emission integrated in one filter. 
Fortunately, this assumption is valid to a good approximation for several existing 
instruments. The passbands of the two Extreme UltraViolet Imagers (EUVI: 



Wuelser et al.\ 2004) on STEREO and those of the Extreme ultraviolet Imaging 



Telescope (FIT: Delaboudiniere et al. 1995 I have for example been designed to 
be identical. 

Equation (J9j) can be reformulated as Equation (|T|) by a simple concatenation 



as shown in Equation (10) 



SOLA: siddon.tex; 18 January 2013; 20:30; p. 9 



N. Barbey et al. 



Finally, multiple-spacecraft tomography in the static case can be formulated 
as single-spacecraft tomography as long as all of the instruments have the same 
spectral bandwidth. The only modification is the multiplication by gains, which 
vary from one instrument to the other. We can then write the estimated map 



using Equation M| as in Equation (111 



esRT,A,B = argmin {JIj/a.b - Pa, se|| 2 + A \\DRe 



(11) 



Using a DEM model, it is possible to extend multiple-spacecraft tomography 
to cases in which the spacecraft have different bandpasses since the spectral 
response is then integrated into the model. This is the approach followed by 



(Frazin, Vasquez, and Kamalabadi, 2009). 



3.3.3. Smooth temporal tomography 



Because the corona is not static, dynamic models are desirable. In this case 
however, it is no longer possible to simplify Equation (|6| as in Equation 
When temporal evolution is present, the recording of data taken at different 
times can be expressed as 



y 



\ Vt 




( ni 



P T e T + n (12) 



This results in a highly underdetermined inverse problem since there are 
T times more unknowns than in Equation ([7]). This underdetermination can 
be mitigated using either priors, such as a temporal smoothness prior, or a 
parameterization of the temporal evolution. Thus a classic, smooth, temporal 
solar-rotational tomography (STSRT) would perform conjugate gradient esti- 



mation using the criterion given in Equation (13 1, where D r and D t are the 



finite-difference operators in space and time. Typically, due to strong underde- 
termination, the hyperparameter X t would be greater than A r , favoring solutions 
with small temporal changes. 



esTSRT = || y - Prer II +A r .||D r e r || +A t ||D t er| 



(13) 



This kind of approach using spatio-temporal regularization has been explored 



before (Zhang, Ghodrati, and Brooks, 2005 Khalsa and Fessler, 20071. 



Following the same approach as in Section 3.3.2 it is possible to generalize 
this expression to the case of multiple spacecraft. 



3.4. EUV Lines 



In the case of EUV lines or EUV bands, the dominant formation process of the 
observed radiation is excitation by collisions between ions and electrons. The 
local emissivity can thus be supposed to be isotropic in which case the quantity 
inverted in Sections |3.3.1| to |3.3.3| is directly the local emissivity of the plasma 
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Figure 1. Geometry of Thomson scattering in the corona. M is the location of the scattering 
electron, O is Sun center, S is an emission point on the solar surface, D is the observer, and 
T is the tangent to the solar surface, p is the impact parameter of (DM), r is the distance of 
the scattering point to Sun center, Rq is the solar radius. Q is the angle between the line of 
sight and the tangent to the solar surface passing through the scattering point. 



summed over the spectral response of the instrument. Resonant scattering may 
contribute significantly to EUV bands such as the 17.1 nm and 19.5 nm bands 
used in, e.g., TRACE, EIT, EUVI, and AIA (ISchrijver and McMullen, 2000b. If 



this is confirmed, then the local emissivity is not isotropic and one needs to apply 
a correction factor to the inverted quantities to deduce plasma emissivities. 



3.5. White Light: Thomson Scattering 



In the case of white-light detectors, the measured intensity is largely dominated 
by Thomson scattering of the photospheric radiation by free coronal electrons. 
Figure [T shows the geometry of Thomson scattering in the corona. Following 



Billings (19661, the equations of Thomson scattering in the corona are 



It = 
L = 



7TCT To 

2 r 2 

KSL ISLr, £ 

2 r 2 "" e r ; 



[(1 - u)C 3 + uC 4 ] 
[(l-u)Ci+«C a ] 



(14) 



C 2 
C 3 

c 4 



cos 57 sin Jl 



1 — 3 Sin 2 SI — £°s_£lf1 _L Q oir, 2 o\ 1,-, ( l+sinjp 



sin Q 



(1 + 3 sin 2 fi) In I 



cos n 



cos 



(15) 



where u accounts for the center-to-limb variation and is a function of wavelength, 
7 t and I r are intensities in the radial and transverse directions, n e is the electron 
density, r is the distance of the scattering point to the center of the Sun, p is the 
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impact parameter of the line of sight, Q is the angle between the line of sight 
and the tangent to the solar surface passing through the scattering point, a is 
the Thomson-scattering cross-section and Iq is the incident intensity. 

The important point in these equations is that the intensity is a linear func- 
tion of the electron density. It is thus possible to directly estimate the electron 
density using solar tomography on white-light data. Note also that the Thomson- 
scattering equations can be separated into coefficients that depend on the po- 
sition [M] in the corona (through r and f2) and the line of sight (through the 
impact parameter [p\). These coefficients are given in Equation (16 1. 



Ok 

m : 



Pi 

(l-«)Ci(rj)+C 2 (r.,-) 
(l- M )C 3 (r 3 )+C 4 (r J ) 



(16) 



Measurements are generally decomposed into a polarized brightness [pB] com- 
ponent and a total brightness [B] component. Equation (17 1 gives the equations 
for those quantities. 



pB k = 



B,. 



! J/o|-^[(l-^)C 1 (r,)+ U C 2 (r J )] 
+2^ [(1 ~ u)C 3 ( rj ) + uC 4 { rj )} \ n eJ 



(17) 



From Equations ( 16 ) and ( 17 1 one can build linear direct models for pB and B 



images as in Equation ( 18 1 , where n e is the discretized electron-density map and 



i p B and %b are the pB and B images respectively. O, M, and M' are diagonal 
matrices with o^, rrij, and m'j on their diagonals. 



IB 



T p Bl1 e + fl 

OPMn, + n 

T B n e + n 

P 
P 



O I 



M 

M' 



(18) 



Note that almost twice as much computation is required for simulation or 
inversion of total brightness data than for polarized brightness data. A recent 



application of this model can be found in Frazin et al. (20101 



4. Performance 

We performed tests on a set of 64 images of 256 x 256 pixels, and a reconstruction 
cube of 128 x 128 x 128 voxels. We always use these parameters unless specified 
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Table 1. Performances of the projector as a function of threads. 



Cores 


1 


2 


3 


4 


5 


6 


7 


8 


Projection time [s] 


97.5 


48.6 


36.3 


26.4 


21.2 


18.2 


16.8 


13.9 


back-projection time [s] 


145.3 


76.8 


56.5 


41.0 


37.0 


32.7 


30.1 


23.9 



otherwise. For these tests there is no obstacle, meaning that the ray tracing is 
not interrupted as it would have been with the use of a model with an opaque 
photosphere. There is no mask applied to the data or the map. Tests have been 
done on a PC with two Quad-Core AMD Opteron™ Processor 2380 and 32 
Gigabytes of RAM. 

In a more realistic use of TomograPy for solar tomography, the projections 
would be faster than presented here since masked pixels are not projected and 
LOS integration is stopped when the ray reaches the photosphere, reducing the 
number of computations. 

Performance as a function of the number of threads used by OpenMP are 
shown in Table [T] It shows that the time to compute a projection and a back- 
projection is not exactly linear with the number of threads. It takes seven times 
less time to compute a projection with eight threads than with one thread and 6.1 
times less to compute a back-projection. Back-projection does not scale as well 
as projection (the speed-up with multiple cores is better with the projection). 
This is due to the fact that extra care must be taken while updating the map of 
voxel values as opposed to the detector values in order to avoid race conditions. 
Race conditions are situation in which the outcome of a computation varies 
unexpectedly due to the timing of events in different threads. In TomograPy, 
race conditions occur mainly when two LOS need to update the same voxel at 
the same time. In the case of the projection, voxel values are only read, so this 
not an issue. Race conditions are avoided with the OpenMP atomic pragma 
directive which instruct each thread to update voxels sequentially, resulting in a 
slow-down of computations, but only for this part of the algorithm. Note however 
that the projections and back-projections scale better when the ratio between 
the number of data samples and the number of voxels decreases. Indeed, in this 
case, the number of lines of sight intersecting a single voxel decreases so that it 
is less probable that several threads try to update the same voxel at the same 
time. 

Performance as a function of the image format is shown in Table [2] As 
expected, the projection and back-projection duration scales linearly with the 
image size. 

Performance as a function of the reconstruction cube size is shown in Table [3] 
We can see that the projection and back-projection duration scales linearly with 
the cube root of the number of voxels. This is expected as the number of voxels 
along a single LOS is roughly proportional to the number of voxels along one 
axis of the map and the number of operations is proportional with the number 
of intersections between LOS and voxels. 
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Table 2. Performance of the projector as a function of image shape. 



Image size 


128 X 128 


256 x 256 


512 x 512 


1024 x 1024 


Projection time [s] 


0.95 


3.64 


14.53 


58.75 


back-projection time [s] 


1.61 


6.17 


24.69 


99.54 



Table 3. Performance of the projector as a function of the size 
of the reconstruction cube. 



Cube size 


128 3 


256 3 


512 s 


1024 3 


Projection time [s] 


13.84 


44.56 


195.69 


587.56 


back-projection time [s] 


24.71 


56.64 


214.55 


610.75 



5. Examples 

5.1. Static Reconstruction using STEREO/EUVI A and B 

We performed a conjugate gradient inversion using data from both EUVI A and 
B. To avoid issues due to differences in filters, we rescaled the EUVI B data to 
EUVI A levels by dividing by the empirically deduced values provided in Table 

m 

This operation is justified since the passbands of the twin instruments were 
designed and were measured to be nearly identical. To first order, the difference 
in spectral response between EUVI A and B is a scaling factor. The image pairs 
where chosen to be simultaneous with a small lossy-compression factor. Because 
of the photosphere opacity, we need a full rotation (four weeks) to have data 
isotropic coverage of all parts of the reconstruction cube. If we focus only on 
the poles, this is no longer an issue and we can use only two weeks of data. We 
can also use the diversity of points of view provided by STEREO to reduce this 
duration. At the time of observation, the STEREO spacecraft were separated 
by 86°, reducing the acquisition time needed for a complete coverage of the 
corona to three weeks instead of a full solar rotation. For polar regions, the 
required acquisition time is further reduced to 11 days. We chose pairs of images 
regularly spaced from 1 to 15 December 2008, with four pairs of images per day 
and per observatory, resulting in 118 images. The estimated 3D map is a cube 

Table 4. Ratio of sensitivity between 



EUVI A and 


EUVI 


B as 


a function of 


wavelength. 










bandwidth 


171 


195 


284 


304 


B/A ratio 


0.90 


0.97 


0.95 


1.05 
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DN/m 

Figure 2. Reconstruction of the 17.1 nm local emissivity at 1.05 solar radii using 
STEREO/EUVI A and B data and assuming a static corona. Input data are 59 pairs of 
images equally spaced in time from 1 to 15 December 2008. The STEREO separation angle 
was 86°. 



of 256 x 256 x 256 voxels centered on the Sun with a width of three solar radii 
along each axis. In order to save computation time while remaining consistent 
with the resolution of the reconstruction cube, the images were binned 2x2. 
The Carrington rotation rate is assumed. 

Figure ^ shows the local emissivity in the reconstructed cube at a constant 
altitude of 1.05 solar radii. An equi-rectangular projection is used. 

As can be seen, part of the projection around 300° of longitude is smoother 
than the remainder of the projection. This is due to the fact that we used two 
weeks of data instead of the three weeks required to sample all voxels isotropically 
with lines of sight. To obtain a better estimate at these locations would have 
required the use of one more week's data, which would have in turned worsened 
artifacts due to temporal evolution. Conversely, the polar regions are slightly 
oversampled because, considering the separation of the STEREO spacecraft, the 
minimum integration time was 11 days. Figure [4] shows a gnomonic projection 
of the north polar region at 1.05 solar radii. The coronal hole is clearly visible 
as the darker central area. Structures in the hole are arranged according to a 
network pattern. One can identify bright nodules that could be attributed to 
the classical beam plumes but also elongated structures that could correspond 



to curtain plumes. This would confirm the proposition by Gabriel et al. (2009) 
that both types of plumes coexist. 

As can be seen, some of the voxels have negative values. This is usually 
explained as resulting from temporal evolution. Indeed, if temporal evolution 
has occurred during the acquisition of the data, it cannot be correctly modeled 
with the static assumption made in this inversion. Negative values are thus 
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required to account for a variation of intensity in the data unexplainable by a 
simple change of viewpoint. This could also be explained by niismodeling of the 
measurement process, noise or even bias in the data (which could be due to some 
instrumental artifcat) . 

To show that convergence has indeed been reached with the conjugate gradient 
algorithm we present the criterion as a function of iteration number in Figure [3] 

5.2. Smooth Temporal Evolution with STEREO/EUVI A and B 



The second example is a 3D map estimation using the smooth temporal evolution 
model. In other words, we minimize Equation (13 I using the hyper-parameters 
stated in Table [5] We used the same set of EUVI images used for the static 
estimation. For each pair of EUVI images there is a corresponding instantaneous 
map of dimensions 128 x 128 x 128, resulting in the estimation of approximately 
132 millions of parameters. The estimation took less than eight hours. 

Figure [5] shows gnomonic projections at 1.05 solar radii of the estimated map 
at different instants separated by 40 hours. It is interesting to compare these 
images with the static reconstruction of the same area shown in Figure [4] Here, 
the disappearance of an elongated structure at the south edge of the coronal 
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Table 5. Hyper-parameters used for the 
inversion of the different models. 



model 






A z 


At 


SRT 


le" 1 


le- 1 


le- 1 




Thomson 


le" 1 


le- 1 


le- 1 




STSRT 


le" 1 


le- 1 


le- 1 


le 2 



hole is very clear, and one can also follow the appearance and disappearance of 
beam plumes. 

The static reconstruction is sharper than the temporally evolving one. This 
could be due to slower convergence in the smooth temporal-rotational tomog- 
raphy model. Indeed, since the temporal prior is much higher than the spatial 
prior, numerous very small steps in the conjugate-gradient algorithm could be 
required to reach the minimum. This could be missed by our stopping criterion 
on the norm of the gradient and even on convergence diagonstic such as the one 
in Figure [3] This could be solved through the use of preconditioning but has not 
been tried for now. However, one can clearly identify the sanies structures in both 
reconstructions. One can picture the static reconstruction as a kind of average 
over time, although this is not strictly true as temporal effects and changes of 
viewpoints can have the same kind of effects on data. 



5.3. Thomson scattering with COR1 A and B data 



We estimated the coronal electron density using COR1 A and B data acquired 



during February 2008 as done by Kramar et al. (20091. Since the inversion 



codes are different, the comparison gives an estimate of the robustness of both 
techniques. Our results are shown in Figure [6] and can be compared to Figure 2 
of Kramar et al. ( 2009 1 . The comparison shows that very large scale structures 
are very similar in both maps, but fainter and smaller scales structures differ. 
This can be explained by the use of different prior models and hyper-parameters. 
Data can also differ in the way that they are prepared before the tomographic 
inversion. Note that we used a smoothness prior increasing with height for this 
reconstruction. 



6. Conclusions 



We developed a tested, fast, and flexible program to perform rotational tomog- 
raphy of the solar corona. Its respect of WCS standards allows its use with 
virtually any data set available. 

We demonstrated how this software can be used to perform multi-spacecraft 
estimations of emission maps in the corona using STEREO/EUVI data and 
STEREO/COR1 data. Estimations were performed using static, temporally- 
evolving, and Thomson-scattering models. 
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Figure 5. Polar view of the emissivity at 1.05 solar radii from a temporally evolving 3D map 
estimation using STEREO/EUVI A and B data in the 17.1 nm passband. North pole gnomonic 
projection with a resolution of 1° in square is used. There is 40 hours between each projection. 
These images can be compared with the static reconstruction of Figure [3] 
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Longitude (degrees) 

Figure 6. Equi-rectangular projection of a static 3D map estimation using pB COR1 A and B 
data during February 2008 at 1.6 Rq (a), 2.0 Rq (b), and 2.4 Rq (c). Pixels of the projections 
are 1° square. 



This new software suite can naturally be used to perform tomographic inver- 
sion of the corona using SDO/AIA images. TomograPy will allow tomographic 
estimation using AIA data at full resolution, providing unprecedented resolution 
as well as additional spectral information. Supplemented with spectral-inversion 
methods, this will allow for electron-density and temperature estimates in large 
data sets, close to the photosphere, at very high resolution. 

Because of its modular architecture, TomograPy also provides a convenient 
way to test different estimation algorithms saving the need to rewrite anything 
other than the algorithm. This allows easy comparison of the performances of 
different algorithms for the solar tomography application. 
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